function [hmF2] = hMF2(foE, foF2, M3000)
%
% F2-Layer Maximum Density Height
%
%DESCRIPTION:
%This function computes the F2-Layer maximum density height.
%
%PROTOTYPE:
% [hmF2] = hMF2(foE, foF2, M3000F2)
%
%--------------------------------------------------------------------------
% INPUTS:
%   foE        [1x1]       E-Layer Crit. Frequency   [MHz]
%   foF2       [1x1]       F2-Layer Crit. Frequency  [MHz]
%   M3000      [1x1]       M(3000)F2                 [-]  
%--------------------------------------------------------------------------
% OUTPUTS:
%   hmF2       [1x1]       F2-Lay. Max. Dens. Height [km]
%--------------------------------------------------------------------------
%
%NOTES:
% (none)
%
%CALLED FUNCTIONS:
% (none)
%
%UPDATES:
% (none)
%
%REFERENCES:
% [1] "Ionospheric Correction Algorithm for Galileo Single-Frequency Users"
%      - European GNSS (Galileo) Open Service
%
%AUTHOR(s):
%Luigi De Maria, Matteo D'Addazio, 2022
%

%% Main Code

%Parameters
rho = (foF2/foE * exp(20*(foF2/foE - 1.75)) + 1.75) /...
                 (exp(20*(foF2/foE - 1.75)) + 1);
%M
M = M3000;
if     foE < 1e-30
    DM = -0.012;
elseif foE >= 1e-30
    DM = (0.253) / (rho - 1.215) - 0.012;
end

%F2-Layer Maximum Density Height
hmF2 = (1490*M*sqrt((0.0196*M^2 + 1) / (1.2967*M^2 - 1))) / (M + DM) - 176;

end